##%%%%%%%%%%%%%%% MERGE.R %%%%%%%%%%%%%%%%%%%%%%%%%%%%
##%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
## In this script, we merge the files  
## generated by the script "InOutRsq.R"

library(Matrix)
library(MASS) ## For ginv() - Moore penrose inverse

##%%%%%%%% Set Working Directory %%%%%%%%%%%%
WorkDir="D:/RWORK/Repository Files"


setwd(WorkDir)

##%%%%%%%%%% Read in function files %%%%%%%%%

source("FUNC.R")

##%%%%%%%%%%%% Load and merge rsqA data %%%%%%%%%%%%%%%%%%
##%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
nameSA7=paste("A7",c("t1N5","t2N5","t1N2","t2N2","t1N1","t2N1","t1N.5","t2N.5"),sep="")
fileSA7=paste("rsq",nameSA7,sep="")
nameSA7S=paste("A7S",c("t1N5","t2N5","t1N2","t2N2","t1N1","t2N1","t1N.5","t2N.5"),sep="")
fileSA7S=paste("rsq",nameSA7S,sep="")
nameSA4S=paste("A4S",c("t1N5","t2N5","t1N2","t2N2","t1N1","t2N1","t1N.5","t2N.5"),sep="")
fileSA4S=paste("rsq",nameSA4S,sep="")
nameSA4=paste("A4",c("t1N5","t2N5","t1N2","t2N2","t1N1","t2N1","t1N.5","t2N.5"),sep="")
fileSA4=paste("rsq",nameSA4,sep="")
nameSA3=paste("A3",c("t1N5","t2N5","t1N2","t2N2","t1N1","t2N1","t1N.5","t2N.5"),sep="")
fileSA3=paste("rsq",nameSA3,sep="")
nameSA2=paste("A2",c("t1N5","t2N5","t1N2","t2N2","t1N1","t2N1","t1N.5","t2N.5"),sep="")
fileSA2=paste("rsq",nameSA2,sep="")


## IN the full case don't forget fileSA2, nameSA2!!
fileS=c(fileSA7,fileSA7S,fileSA4S,fileSA4,fileSA3,fileSA2)
nameS=c(nameSA7,nameSA7S,nameSA4S,nameSA4,nameSA3,nameSA2)

## Load files
for(i in 1:length(fileS)){ load(fileS[i]) }

## Merge files
foo=get(fileS[1])
rownames(foo)=paste(nameS[1],rownames(foo))
rsqAA=foo

for(i in 2:length(fileS)){
	foo=get(fileS[i])
	rownames(foo)=paste(nameS[i],rownames(foo))
	rsqAA=rbind(rsqAA,foo)
	}
	
##%%%%%%%% Get rid of the 600 or so with 0.1 
## large deviations between "t value" and "tstatF"
## even though it changes next to nothing in the 
## analysis.

jj=abs(abs(rsqAA[,"t value"])-rsqAA[,"tstatF"])<0.1

## > sum(jj)/length(jj)
## [1] 0.9946572

rsqAA=rsqAA[jj,]


##%%%%%%%%%%% Indentify sims models and inputs %%%%%%%%%%%%
##%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
rNames=rownames(rsqAA)
simNames=substring(rNames,1,8)

fuN=function(strG){
	ouT=strG
	lasT=substring(ouT,nchar(ouT),nchar(ouT))
	ii=lasT=="m"
	ouT[ii]=substring(ouT[ii],1,nchar(ouT[ii])-2)
	ii=lasT==" "
	ouT[ii]=substring(ouT[ii],1,nchar(ouT[ii])-1)
	ouT
	}

simNames=fuN(simNames)
simNamesU=unique(simNames)

modNames=substring(rNames,nchar(simNames)+2,nchar(simNames)+7)
inpNames=substring(rNames,nchar(simNames)+9,nchar(rNames))


##%%%%%%%%%%%%% Compute median/std/mean across models %%%%%%%%%%
##%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
## Get simNamesU ordered by n-obs
ordStrgs=c("N5","N2","N1",".5")
simNamesUO=c(simNamesU[substring(simNamesU,nchar(simNamesU)-1,nchar(simNamesU))==ordStrgs[1]],
		simNamesU[substring(simNamesU,nchar(simNamesU)-1,nchar(simNamesU))==ordStrgs[2]],
		simNamesU[substring(simNamesU,nchar(simNamesU)-1,nchar(simNamesU))==ordStrgs[3]],
		simNamesU[substring(simNamesU,nchar(simNamesU)-1,nchar(simNamesU))==ordStrgs[4]])


##%%%%%% Compute cross model medians for all inputs/stats %%%%%%%%

ii=(simNames==simNamesUO[1])
inpN1=inpNames[ii]
rsqLS1=rsqAA[ii,]

foo=inpStatsFunc(rsqLS1,inpN1,Cmedian)
rownames(foo)=paste(simNamesUO[1],rownames(foo))
rsqInpMedAA=foo

for(j in 2:length(simNamesUO)){
	ii=(simNames==simNamesUO[j])
	inpN1=inpNames[ii]
	rsqLS1=rsqAA[ii,]
	foo=inpStatsFunc(rsqLS1,inpN1,Cmedian)
	rownames(foo)=paste(simNamesUO[j],rownames(foo))
	rsqInpMedAA=rbind(rsqInpMedAA,foo)
	}



##%%%%%% Collect stddev's for all inputs/stats
ii=(simNames==simNamesUO[1])
inpN1=inpNames[ii]
rsqLS1=rsqAA[ii,]

foo=inpStatsFunc(rsqLS1,inpN1,Cstd)
rownames(foo)=paste(simNamesUO[1],rownames(foo))
rsqInpStdAA=foo

for(j in 2:length(simNamesUO)){
	ii=(simNames==simNamesUO[j])
	inpN1=inpNames[ii]
	rsqLS1=rsqAA[ii,]
	foo=inpStats(rsqLS1,inpN1,Cstd)
	rownames(foo)=paste(simNamesUO[j],rownames(foo))
	rsqInpStdAA=rbind(rsqInpStdAA,foo)
	}



##%%%%%%%% Create rsqMedAA all medians matched to inputs %%%%%%%%%%%%%
## rsqMedAA will be the same size as rsqAA, but will only contain
## cross model median's of the relevant sim/input combinations. This 
## will be used when comparing individual model stats with medians
## across models.
simInpNames=paste(simNames,inpNames)
simInpNamesS=rownames(rsqInpMedAA)

rsqMedAA=cbind(matrix(1,nrow(rsqAA),1),rsqAA-rsqAA)
colnames(rsqMedAA)=colnames(rsqInpMedAA)

for(i in 1:length(simInpNamesS)){
	simInpN=simInpNamesS[i]
	jj=(simInpNames==simInpN)
	rsqMedAA[jj,]=matrix(1,sum(jj),1)%*%rsqInpMedAA[simInpN,]
	}


## Do the same thing for standard deviation
rsqStdAA=cbind(matrix(1,nrow(rsqAA),1),rsqAA-rsqAA)
colnames(rsqStdAA)=colnames(rsqInpStdAA)

for(i in 1:length(simInpNamesS)){
	simInpN=simInpNamesS[i]
	jj=simInpNames==simInpN
	rsqStdAA[jj,]=matrix(1,sum(jj),1)%*%rsqInpStdAA[simInpN,]
	}

##%%%%% Add vector of the number of n-pval<=0.1 for each sim to rsqAA %%%%%
rsqAA=cbind(rsqMedAA[,"n-pval<=0.1"],rsqAA)
colnames(rsqAA)[1]=c("n-pval<=0.1")


##%%%%%%%%%%%%%% Save the merged Files %%%%%%%%%%%%%%%%%%%%
# save(rsqAA,file="rsqAA")
# save(rsqInpMedAA,file="rsqInpMedAA")
# save(rsqInpStdAA,file="rsqInpStdAA")
# save(rsqMedAA,file="rsqMedAA")
# save(rsqStdAA,file="rsqStdAA")

